function out = mcmc(SET, data, zlb_t, c, x0, runs, cont_xx)
% 
% Runs the Markov Chain Monte Carlo
%
% Outputs the chains for parameters and durations

if nargin < 6 ; runs = 1 ; end

%% Parameters

nparam = length(x0);
ss     = SET.ss ;

N     = SET.EST.N ; 
df    = SET.EST.df ; 
Tbstar = SET.EST.Tbstar ; % Maximum expected duration at the ZLB is Tbstar;

pop = 1:nparam ;

kk_xx       = 0 ; acc_xx  = 0 ;

%% Posterior

post     = @(x, T_) ll(x, T_, SET, data) + prior(x) ;
post_neg = @(x, T_) -post(x, T_) ;

%% Settings for parameters

D        = sum(zlb_t) ; 
T        = zeros(1,SET.ss) ;

Prop_xx      = zeros(N,length(x0)) ;
Prop_xx(1,:) = x0' ; 
xx_out       = zeros(N,length(x0)) ;
xx_out(1,:)  = Prop_xx(1,:)';

%% If continuing a chain

if nargin>6  
    Prop_xx(1,:) = cont_xx(end,:,runs) ;
    xx_out(1,:)  = cont_xx(end,:,runs) ;
    
    T(zlb_t>0) =  SET.T_f ;
    
    denominator = post_neg(xx_out(1,:), T) ;
end

%% If new chain

if nargin<7
    if SET.end_date>2008.75
        T(zlb_t>0) =  SET.T_f ;
    end

    denominator = post_neg(xx_out(1,:), T) ;
end

%% DUAL SINGLE-MOVE SAMPLER
for j=2:N
%% xx-block

xx_prop     = xx_out(j-1,:);
blocksize   = ceil(nparam*rand);
block       = randsample(pop,blocksize);

xx_prop(1,block) = xx_out(j-1,block) + ...
    (c(block,block)*trnd(df,blocksize,1))'; % multivariate-t with df degrees of freedom

Prop_xx(j,:) = xx_prop(1,:);

numerator  = post_neg(xx_prop, T) ;
ratio2     = exp(-numerator + denominator) ;
alpha2     = min(ratio2,1) ;

if isnan(ratio2) ; alpha2 = 0; end

if rand <= alpha2
   xx_out(j,:) = xx_prop;
   kk_xx       = kk_xx + 1;
   denominator=numerator ;
else
    xx_out(j,:) = xx_out(j-1,:);
end

acc_xx= 100*kk_xx/j;

if mod(j/10,1)==0 
    fprintf(sprintf('Chain %1.0f .... %2.2f per cent done\n   ***** accept rates: x, %2.3f\n',runs,100*j/N,acc_xx)) ;
end

if mod(j/100,1)==0 
    mean(xx_out(1:j,:))
    denominator
end

if mod(j/(N/SET.EST.numbersave),1)==0
    eval(['save chain_', int2str(runs), ';']) ;
end

end
%% Outputs

out.xx_out  = xx_out ;
out.acc_rates = [acc_xx, N] ;